% Script to read and plot 1 hour geomag data from BGS

S = dir('/Users/manojnair/data/obs_mag_data/BGS data/*.nc');

time_array_hour = (datenum(2013,1,1,0,0,30): (1/(24)): datenum(2014,12,31,23,59,30))';

%% plot

f = figure(1);
set(f,'Position',[40          19        1557         937]);
for i = 1:length(S),
    
    ncfname = ['/Users/manojnair/data/obs_mag_data/BGS data/' S(i).name];
    
    
    [x_data, y_data, z_data , X_ID, Y_ID, Z_ID, obj] = read_geomag_netcdf(ncfname, 157800, 17520, 0);
    
    if any(x_data),
    
    if     obj.geospatial_lon > 180,
        obj.geospatial_lon = obj.geospatial_lon - 360;
    end;
    
    icount = 1;
    for decyear = 2013:0.1:2014.5,
        [mag_igrf11(icount,:)] = igrf11magm(0,obj.geospatial_lat,obj.geospatial_lon,decyear);
        [mag_igrf12(icount,:)] = magsynth(0,obj.geospatial_lat,obj.geospatial_lon,decyear);

        icount = icount + 1;
    end;
    
    
    subplot(311);
    
    plot(time_array_hour,x_data);
    hold on;
    plot(datenum(2013:0.1:2014.5,0,0), mag_igrf11(:,1),'r.-','MarkerSize',15);
    plot(datenum(2013:0.1:2014.5,0,0), mag_igrf12(:,1),'b.-','MarkerSize',15);

    legend('OBS','IGRF11','IGRF12');
    set(gca,'FontSize',16);
    
    ylabel('X nT');
    datetick('x',2);
    axis([datenum(2013,1,1), datenum(2015,1,1), -inf inf]);
    title(S(i).name(1:3));
    
    subplot(312);
    
    plot(time_array_hour,y_data);
    hold on;
    plot(datenum(2013:0.1:2014.5,0,0), mag_igrf11(:,2),'r.-','MarkerSize',20);
    plot(datenum(2013:0.1:2014.5,0,0), mag_igrf12(:,2),'b.-','MarkerSize',20);

    
    set(gca,'FontSize',16);
    ylabel('Y nT');
    datetick('x',2);
    axis([datenum(2013,1,1), datenum(2015,1,1), -inf inf])
    subplot(313);
    
    plot(time_array_hour,z_data);
    hold on;
    plot(datenum(2013:0.1:2014.5,0,0), mag_igrf11(:,3),'r.-','MarkerSize',20);
        plot(datenum(2013:0.1:2014.5,0,0), mag_igrf12(:,3),'b.-','MarkerSize',20);

    
    set(gca,'FontSize',16);
    ylabel('Z nT');
    datetick('x',2);
    axis([datenum(2013,1,1), datenum(2015,1,1), -inf inf])
    
    
    set(gcf, 'PaperPositionMode', 'auto');
    %pause;
    saveas(gcf,['/Users/manojnair/projects/igrf12_validation/IGRF12_comp' S(i).name(1:3)],'jpeg');
    clf
    end;
end;


